"""
FF: The Demographic Regulatory Doom Loop

Empirical phases:
1. Data assembly (merge cliff panel + 140-country panel + pension data)
2. Pension endogeneity (OADR → pension spending/coverage)
3. Fiscal-pension feedback (pension spending → debt dynamics)
4. Sovereign risk amplification (pension × demographics × fiscal → spreads/ratings)
5. Doom loop dynamics (local projections: shocks propagate through the loop)
6. Scenario projections (forward-looking doom loop intensity by country)
"""

import sys
sys.path.insert(0, '/mnt/c/demographics_capital_flows/multilateral/140_country/src')

import pandas as pd
import numpy as np
from model import PanelGLS
from pathlib import Path
from scipy import stats

OUTPUT_DIR = Path('/mnt/c/demographics_capital_flows/doom_loop/output/tables')
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

# ══════════════════════════════════════════════════════════════════════════════
# PHASE 1: DATA ASSEMBLY
# ══════════════════════════════════════════════════════════════════════════════
print("=" * 80)
print("PHASE 1: DATA ASSEMBLY")
print("=" * 80)

# Start with cliff panel (31 sovereign-rated countries, 1990-2024)
cliff = pd.read_csv('/mnt/c/demographics_capital_flows/safe_asset_cliff/data/processed/cliff_panel.csv')
cliff = cliff[cliff['year'] <= 2024].copy()

# Merge with 140-country panel for additional variables
full = pd.read_csv('/mnt/c/demographics_capital_flows/multilateral/140_country/data/processed/full_panel.csv')
full = full[full['year'] <= 2024].copy()

# Take pension and additional fiscal vars from full panel
extra_cols = ['iso3', 'year', 'trade_openness']
extra = full[extra_cols].drop_duplicates()

df = cliff.copy()
print(f"Cliff panel: {df.shape[0]} obs, {df['iso3'].nunique()} countries, {df['year'].min()}-{df['year'].max()}")

# Create key derived variables
df['oadr'] = df['old_dep']  # alias
df['log_gdp_pc'] = np.log(df['gdp_pc_ppp'].clip(lower=100))
df['debt_x_oadr'] = df['govt_debt_gdp'] * df['old_dep']
df['pension_x_oadr'] = df['pension_spending_gdp'] * df['old_dep']
df['exp_gap_x_oadr'] = df['exp_rev_gap'] * df['old_dep']
df['safe_x_oadr'] = df['safe_issuer'] * df['old_dep']

# OECD indicator
OECD = {'AUS','AUT','BEL','CAN','CHL','COL','CRI','CZE','DNK','EST','FIN','FRA',
        'DEU','GRC','HUN','ISL','IRL','ISR','ITA','JPN','KOR','LVA','LTU','LUX',
        'MEX','NLD','NZL','NOR','POL','PRT','SVK','SVN','ESP','SWE','CHE','TUR',
        'GBR','USA'}
df['oecd'] = df['iso3'].isin(OECD).astype(float)

# Lagged variables for dynamics
df = df.sort_values(['iso3', 'year'])
for var in ['rating_numeric', 'govt_debt_gdp', 'pension_spending_gdp', 'sovereign_spread',
            'exp_rev_gap', 'old_dep']:
    df[f'{var}_lag'] = df.groupby('iso3')[var].shift(1)
    df[f'{var}_change'] = df[var] - df[f'{var}_lag']

# Rating deterioration indicator
df['rating_deteriorated'] = (df['rating_numeric_change'] < 0).astype(float)

# Debt acceleration
df['debt_acceleration'] = df['govt_debt_gdp_change'] - df.groupby('iso3')['govt_debt_gdp_change'].shift(1)

# Save assembled panel
df.to_csv('/mnt/c/demographics_capital_flows/doom_loop/data/processed/doom_loop_panel.csv', index=False)
print(f"Assembled panel saved: {df.shape[0]} obs, {df['iso3'].nunique()} countries")
print(f"Pension spending coverage: {df['pension_spending_gdp'].notna().sum()} obs")
print(f"Sovereign spread coverage: {df['sovereign_spread'].notna().sum()} obs")
print(f"Rating coverage: {df['rating_numeric'].notna().sum()} obs")


# ══════════════════════════════════════════════════════════════════════════════
# PHASE 2: PENSION ENDOGENEITY — DOES AGING PREDICT PENSION GENEROSITY?
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n{'=' * 80}")
print("PHASE 2: PENSION ENDOGENEITY")
print("=" * 80)

def stars(p):
    if p < 0.01: return '***'
    elif p < 0.05: return '**'
    elif p < 0.10: return '*'
    return ''

def run_gls(dv_name, xvar_names, data, label=""):
    """Run PanelGLS and return results dict."""
    cols = ['iso3', 'year', dv_name] + xvar_names
    sub = data[cols].dropna()
    if len(sub) < 30:
        print(f"  SKIP {label}: only {len(sub)} obs")
        return None
    gls = PanelGLS()
    gls.fit(sub[dv_name].values, sub[xvar_names].values, sub['iso3'].values, sub['year'].values)
    results = {
        'label': label,
        'n': gls.n_obs,
        'r2': gls.r_squared,
        'n_countries': gls.n_countries,
    }
    for j, name in enumerate(xvar_names):
        results[f'{name}_beta'] = gls.beta[j]
        results[f'{name}_se'] = gls.se[j]
        results[f'{name}_p'] = gls.pvalues[j]
    return results

# Test: OADR → pension spending
controls_fiscal = ['log_gdp_pc', 'rgdp_growth', 'inflation']
pension_results = []

# M1: old_dep → pension_spending
r = run_gls('pension_spending_gdp', ['old_dep'] + controls_fiscal, df, "OADR → Pension spending")
if r:
    pension_results.append(r)
    print(f"  old_dep → pension_spending: β={r['old_dep_beta']:.3f} (p={r['old_dep_p']:.4f}{stars(r['old_dep_p'])}), n={r['n']}, R²={r['r2']:.3f}")

# M2: old_dep → health_exp
r = run_gls('health_exp_gdp', ['old_dep'] + controls_fiscal, df, "OADR → Health expenditure")
if r:
    pension_results.append(r)
    print(f"  old_dep → health_exp: β={r['old_dep_beta']:.3f} (p={r['old_dep_p']:.4f}{stars(r['old_dep_p'])}), n={r['n']}, R²={r['r2']:.3f}")

# M3: old_dep → govt_expenditure
r = run_gls('govt_expenditure_gdp', ['old_dep'] + controls_fiscal, df, "OADR → Govt expenditure")
if r:
    pension_results.append(r)
    print(f"  old_dep → govt_exp: β={r['old_dep_beta']:.3f} (p={r['old_dep_p']:.4f}{stars(r['old_dep_p'])}), n={r['n']}, R²={r['r2']:.3f}")

# M4: old_dep → govt_revenue
r = run_gls('govt_revenue_gdp', ['old_dep'] + controls_fiscal, df, "OADR → Govt revenue")
if r:
    pension_results.append(r)
    print(f"  old_dep → govt_rev: β={r['old_dep_beta']:.3f} (p={r['old_dep_p']:.4f}{stars(r['old_dep_p'])}), n={r['n']}, R²={r['r2']:.3f}")

# M5: old_dep → expenditure-revenue gap
r = run_gls('exp_rev_gap', ['old_dep'] + controls_fiscal, df, "OADR → Exp-Rev gap")
if r:
    pension_results.append(r)
    print(f"  old_dep → exp_rev_gap: β={r['old_dep_beta']:.3f} (p={r['old_dep_p']:.4f}{stars(r['old_dep_p'])}), n={r['n']}, R²={r['r2']:.3f}")

# M6: old_dep → debt/GDP
r = run_gls('govt_debt_gdp', ['old_dep'] + controls_fiscal, df, "OADR → Debt/GDP")
if r:
    pension_results.append(r)
    print(f"  old_dep → debt/GDP: β={r['old_dep_beta']:.3f} (p={r['old_dep_p']:.4f}{stars(r['old_dep_p'])}), n={r['n']}, R²={r['r2']:.3f}")

# Write pension endogeneity table
lines = ['# Table 1: Pension Endogeneity — Does Aging Drive Fiscal Expansion?\n']
lines.append('PanelGLS with AR(1). Sample: sovereign-rated issuers, 1990-2024.')
lines.append('Controls: log GDP/pc, GDP growth, inflation.\n')
lines.append('| Dependent Variable | OADR coef | SE | p-value | N | R² |')
lines.append('|:---|---:|---:|---:|---:|---:|')
for r in pension_results:
    dv = r['label'].split('→')[1].strip()
    lines.append(f"| {dv} | {r['old_dep_beta']:.3f} | {r['old_dep_se']:.3f} | {r['old_dep_p']:.4f}{stars(r['old_dep_p'])} | {r['n']} | {r['r2']:.3f} |")
lines.append('')
(OUTPUT_DIR / 'table1_pension_endogeneity.md').write_text('\n'.join(lines))


# ══════════════════════════════════════════════════════════════════════════════
# PHASE 3: FISCAL-PENSION FEEDBACK — DOES PENSION SPENDING WORSEN DEBT?
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n{'=' * 80}")
print("PHASE 3: FISCAL-PENSION FEEDBACK")
print("=" * 80)

feedback_results = []

# M1: pension → debt change
r = run_gls('govt_debt_gdp_change', ['pension_spending_gdp', 'rgdp_growth', 'r_minus_g', 'govt_debt_gdp_lag'], df, "Pension → ΔDebt")
if r:
    feedback_results.append(r)
    print(f"  pension → ΔDebt: β={r['pension_spending_gdp_beta']:.3f} (p={r['pension_spending_gdp_p']:.4f}{stars(r['pension_spending_gdp_p'])})")

# M2: pension × oadr → debt change (doom loop interaction)
df['pension_x_old_dep'] = df['pension_spending_gdp'] * df['old_dep']
r = run_gls('govt_debt_gdp_change', ['pension_spending_gdp', 'old_dep', 'pension_x_old_dep', 'rgdp_growth', 'r_minus_g', 'govt_debt_gdp_lag'], df, "Pension×OADR → ΔDebt")
if r:
    feedback_results.append(r)
    print(f"  pension×OADR → ΔDebt: interaction β={r['pension_x_old_dep_beta']:.3f} (p={r['pension_x_old_dep_p']:.4f}{stars(r['pension_x_old_dep_p'])})")

# M3: exp_rev_gap → debt change
r = run_gls('govt_debt_gdp_change', ['exp_rev_gap', 'rgdp_growth', 'r_minus_g', 'govt_debt_gdp_lag'], df, "Exp-Rev gap → ΔDebt")
if r:
    feedback_results.append(r)
    print(f"  exp_rev_gap → ΔDebt: β={r['exp_rev_gap_beta']:.3f} (p={r['exp_rev_gap_p']:.4f}{stars(r['exp_rev_gap_p'])})")

# M4: exp_rev_gap × oadr → debt change
r = run_gls('govt_debt_gdp_change', ['exp_rev_gap', 'old_dep', 'exp_gap_x_oadr', 'rgdp_growth', 'r_minus_g', 'govt_debt_gdp_lag'], df, "ExpGap×OADR → ΔDebt")
if r:
    feedback_results.append(r)
    print(f"  exp_gap×OADR → ΔDebt: interaction β={r['exp_gap_x_oadr_beta']:.3f} (p={r['exp_gap_x_oadr_p']:.4f}{stars(r['exp_gap_x_oadr_p'])})")

# Write feedback table
lines = ['# Table 2: Fiscal-Pension Feedback\n']
lines.append('Does pension spending accelerate debt accumulation, especially in aging economies?\n')
lines.append('| Model | Key Variable | coef | SE | p-value | N | R² |')
lines.append('|:---|:---|---:|---:|---:|---:|---:|')
for r in feedback_results:
    label = r['label']
    # Find the key variable (first non-control)
    key_var = label.split('→')[0].strip().split('×')[0].strip().lower().replace(' ', '_')
    for vname in ['pension_x_old_dep', 'exp_gap_x_oadr', 'pension_spending_gdp', 'exp_rev_gap']:
        if f'{vname}_beta' in r:
            lines.append(f"| {label} | {vname} | {r[f'{vname}_beta']:.3f} | {r[f'{vname}_se']:.3f} | {r[f'{vname}_p']:.4f}{stars(r[f'{vname}_p'])} | {r['n']} | {r['r2']:.3f} |")
            break
lines.append('')
(OUTPUT_DIR / 'table2_fiscal_pension_feedback.md').write_text('\n'.join(lines))


# ══════════════════════════════════════════════════════════════════════════════
# PHASE 4: SOVEREIGN RISK AMPLIFICATION — THE DOOM LOOP TEST
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n{'=' * 80}")
print("PHASE 4: SOVEREIGN RISK AMPLIFICATION")
print("=" * 80)

doom_results = []

# M1: Baseline rating model (replicating safe asset cliff)
base_controls = ['govt_debt_gdp', 'rgdp_growth', 'inflation', 'log_gdp_pc']
r = run_gls('rating_numeric', ['old_dep'] + base_controls, df, "Baseline: OADR → rating")
if r:
    doom_results.append(r)
    print(f"  Baseline OADR → rating: β={r['old_dep_beta']:.3f} (p={r['old_dep_p']:.4f}{stars(r['old_dep_p'])}), R²={r['r2']:.3f}")

# M2: Add pension spending
r = run_gls('rating_numeric', ['old_dep', 'pension_spending_gdp'] + base_controls, df, "OADR + Pension → rating")
if r:
    doom_results.append(r)
    print(f"  + pension: OADR β={r['old_dep_beta']:.3f} (p={r['old_dep_p']:.4f}), pension β={r['pension_spending_gdp_beta']:.3f} (p={r['pension_spending_gdp_p']:.4f})")

# M3: Add exp-rev gap
r = run_gls('rating_numeric', ['old_dep', 'exp_rev_gap'] + base_controls, df, "OADR + ExpGap → rating")
if r:
    doom_results.append(r)
    print(f"  + exp_gap: OADR β={r['old_dep_beta']:.3f} (p={r['old_dep_p']:.4f}), gap β={r['exp_rev_gap_beta']:.3f} (p={r['exp_rev_gap_p']:.4f})")

# M4: Doom loop triple interaction: OADR × debt × exp_gap
df['doom_triple'] = df['old_dep'] * df['govt_debt_gdp'] * df['exp_rev_gap'] / 1000  # scale
r = run_gls('rating_numeric', ['old_dep', 'govt_debt_gdp', 'exp_rev_gap', 'doom_triple', 'rgdp_growth', 'inflation', 'log_gdp_pc'], df,
            "Doom loop: OADR×Debt×ExpGap → rating")
if r:
    doom_results.append(r)
    print(f"  Doom triple → rating: β={r['doom_triple_beta']:.4f} (p={r['doom_triple_p']:.4f}{stars(r['doom_triple_p'])})")

# M5: OADR × debt interaction on rating
r = run_gls('rating_numeric', ['old_dep', 'govt_debt_gdp', 'debt_x_oadr', 'exp_rev_gap', 'rgdp_growth', 'inflation', 'log_gdp_pc'], df,
            "OADR×Debt → rating")
if r:
    doom_results.append(r)
    print(f"  OADR×Debt → rating: β={r['debt_x_oadr_beta']:.4f} (p={r['debt_x_oadr_p']:.4f}{stars(r['debt_x_oadr_p'])})")

# M6: Rating on sovereign spread
r = run_gls('sovereign_spread', ['old_dep', 'govt_debt_gdp', 'debt_x_oadr', 'exp_rev_gap', 'rgdp_growth', 'log_gdp_pc'], df,
            "OADR×Debt → spread")
if r:
    doom_results.append(r)
    print(f"  OADR×Debt → spread: β={r['debt_x_oadr_beta']:.4f} (p={r['debt_x_oadr_p']:.4f}{stars(r['debt_x_oadr_p'])})")

# M7: Pension spending on sovereign spread
r = run_gls('sovereign_spread', ['old_dep', 'pension_spending_gdp', 'pension_x_old_dep', 'govt_debt_gdp', 'rgdp_growth', 'log_gdp_pc'], df,
            "Pension×OADR → spread")
if r:
    doom_results.append(r)
    print(f"  Pension×OADR → spread: β={r.get('pension_x_old_dep_beta', 'N/A')}")

# Write doom loop table
lines = ['# Table 3: Sovereign Risk Amplification — The Doom Loop\n']
lines.append('Does the combination of aging + pension spending + debt accelerate sovereign risk?\n')
lines.append('| Model | N | R² | Key Finding |')
lines.append('|:---|---:|---:|:---|')
for r in doom_results:
    # Summarize key finding
    key_parts = []
    for vname in ['old_dep', 'pension_spending_gdp', 'exp_rev_gap', 'doom_triple', 'debt_x_oadr', 'pension_x_old_dep']:
        if f'{vname}_beta' in r:
            key_parts.append(f"{vname}={r[f'{vname}_beta']:.3f} (p={r[f'{vname}_p']:.3f})")
    lines.append(f"| {r['label']} | {r['n']} | {r['r2']:.3f} | {'; '.join(key_parts[:2])} |")
lines.append('')
(OUTPUT_DIR / 'table3_doom_loop_amplification.md').write_text('\n'.join(lines))


# ══════════════════════════════════════════════════════════════════════════════
# PHASE 5: LOCAL PROJECTIONS — DOOM LOOP DYNAMICS
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n{'=' * 80}")
print("PHASE 5: LOCAL PROJECTIONS — DOOM LOOP DYNAMICS")
print("=" * 80)

# Test: does a pension spending shock propagate to debt and then to ratings?
# Simple local projection: y_{t+h} - y_t = α + β * shock_t + controls + ε

lp_results = []

for h in [1, 2, 3, 5]:
    # Create forward-looking change in rating
    df[f'rating_change_{h}'] = df.groupby('iso3')['rating_numeric'].shift(-h) - df['rating_numeric']
    df[f'debt_change_{h}'] = df.groupby('iso3')['govt_debt_gdp'].shift(-h) - df['govt_debt_gdp']
    df[f'spread_change_{h}'] = df.groupby('iso3')['sovereign_spread'].shift(-h) - df['sovereign_spread']

    # LP1: pension shock → rating change at horizon h
    r = run_gls(f'rating_change_{h}',
                ['pension_spending_gdp', 'old_dep', 'pension_x_old_dep', 'govt_debt_gdp', 'rgdp_growth', 'log_gdp_pc'],
                df, f"Pension×OADR → Δrating(h={h})")
    if r and 'pension_x_old_dep_beta' in r:
        lp_results.append({
            'h': h, 'dv': 'rating',
            'beta': r['pension_x_old_dep_beta'],
            'se': r['pension_x_old_dep_se'],
            'p': r['pension_x_old_dep_p'],
            'n': r['n']
        })
        print(f"  h={h}: pension×OADR → Δrating: β={r['pension_x_old_dep_beta']:.4f} (p={r['pension_x_old_dep_p']:.4f}{stars(r['pension_x_old_dep_p'])})")

    # LP2: exp_gap shock → debt change at horizon h
    r = run_gls(f'debt_change_{h}',
                ['exp_rev_gap', 'old_dep', 'exp_gap_x_oadr', 'govt_debt_gdp', 'rgdp_growth', 'log_gdp_pc'],
                df, f"ExpGap×OADR → ΔDebt(h={h})")
    if r and 'exp_gap_x_oadr_beta' in r:
        lp_results.append({
            'h': h, 'dv': 'debt',
            'beta': r['exp_gap_x_oadr_beta'],
            'se': r['exp_gap_x_oadr_se'],
            'p': r['exp_gap_x_oadr_p'],
            'n': r['n']
        })
        print(f"  h={h}: exp_gap×OADR → ΔDebt: β={r['exp_gap_x_oadr_beta']:.4f} (p={r['exp_gap_x_oadr_p']:.4f}{stars(r['exp_gap_x_oadr_p'])})")

# Write LP table
lines = ['# Table 4: Local Projections — Doom Loop Dynamics\n']
lines.append('Does the fiscal-demographic interaction propagate forward?\n')
lines.append('| Horizon | DV | Interaction β | SE | p-value | N |')
lines.append('|---:|:---|---:|---:|---:|---:|')
for r in lp_results:
    lines.append(f"| h={r['h']} | Δ{r['dv']} | {r['beta']:.4f} | {r['se']:.4f} | {r['p']:.4f}{stars(r['p'])} | {r['n']} |")
lines.append('')
(OUTPUT_DIR / 'table4_local_projections.md').write_text('\n'.join(lines))


# ══════════════════════════════════════════════════════════════════════════════
# PHASE 6: COUNTRY-LEVEL DOOM LOOP INTENSITY SCORES
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n{'=' * 80}")
print("PHASE 6: COUNTRY-LEVEL DOOM LOOP INTENSITY")
print("=" * 80)

# Create a composite doom loop score for each country in 2024 or latest year
latest = df.groupby('iso3').apply(lambda g: g.dropna(subset=['old_dep','govt_debt_gdp','exp_rev_gap']).iloc[-1] if len(g.dropna(subset=['old_dep','govt_debt_gdp','exp_rev_gap'])) > 0 else None)
latest = latest.dropna(how='all').reset_index(drop=True)

if len(latest) > 0:
    # Normalize components to 0-1
    for var in ['old_dep', 'govt_debt_gdp', 'exp_rev_gap', 'pension_spending_gdp']:
        if var in latest.columns and latest[var].notna().sum() > 0:
            vmin, vmax = latest[var].min(), latest[var].max()
            if vmax > vmin:
                latest[f'{var}_norm'] = (latest[var] - vmin) / (vmax - vmin)
            else:
                latest[f'{var}_norm'] = 0.5

    # Doom loop score: weighted average of demographic pressure, debt, expenditure gap
    weights = {'old_dep_norm': 0.30, 'govt_debt_gdp_norm': 0.30, 'exp_rev_gap_norm': 0.25}
    if 'pension_spending_gdp_norm' in latest.columns:
        weights['pension_spending_gdp_norm'] = 0.15
    else:
        # redistribute
        weights = {'old_dep_norm': 0.35, 'govt_debt_gdp_norm': 0.35, 'exp_rev_gap_norm': 0.30}

    latest['doom_score'] = 0
    for var, w in weights.items():
        if var in latest.columns:
            latest['doom_score'] += w * latest[var].fillna(0)

    # Also get projected OADR from UN projections (use oadr_plus20 as proxy for forward pressure)
    if 'oadr_plus20' in latest.columns:
        latest['forward_oadr'] = latest['oadr_plus20']

    # Sort by doom score
    result_cols = ['iso3', 'year', 'old_dep', 'govt_debt_gdp', 'exp_rev_gap',
                   'pension_spending_gdp', 'rating_numeric', 'safe_issuer', 'doom_score']
    result_cols = [c for c in result_cols if c in latest.columns]
    result = latest[result_cols].sort_values('doom_score', ascending=False)

    print("\nDoom Loop Intensity Rankings (latest data):")
    print("-" * 80)
    for _, row in result.head(15).iterrows():
        pension_str = f"{row.get('pension_spending_gdp', np.nan):.1f}" if pd.notna(row.get('pension_spending_gdp')) else 'N/A'
        rating = int(row['rating_numeric']) if pd.notna(row.get('rating_numeric')) else 'N/A'
        safe = 'Yes' if row.get('safe_issuer', 0) == 1 else 'No'
        print(f"  {row['iso3']:5s}  OADR={row['old_dep']:.1f}%  Debt={row['govt_debt_gdp']:.0f}%  Gap={row['exp_rev_gap']:.1f}pp  "
              f"Pension={pension_str}%  Rating={rating}  Safe={safe}  Score={row['doom_score']:.3f}")

    # Write country table
    lines = ['# Table 5: Doom Loop Intensity by Country\n']
    lines.append('Composite score based on demographic pressure (OADR), debt/GDP, expenditure-revenue gap,')
    lines.append('and pension spending. Higher = more vulnerable to doom loop dynamics.\n')
    lines.append('| Country | OADR | Debt/GDP | Exp-Rev Gap | Pension/GDP | Rating | Safe | Doom Score |')
    lines.append('|:---|---:|---:|---:|---:|---:|:---|---:|')
    for _, row in result.iterrows():
        pension_str = f"{row.get('pension_spending_gdp', np.nan):.1f}" if pd.notna(row.get('pension_spending_gdp')) else '—'
        rating = f"{int(row['rating_numeric'])}" if pd.notna(row.get('rating_numeric')) else '—'
        safe = 'Yes' if row.get('safe_issuer', 0) == 1 else 'No'
        lines.append(f"| {row['iso3']} | {row['old_dep']:.1f} | {row['govt_debt_gdp']:.0f} | {row['exp_rev_gap']:.1f} | "
                     f"{pension_str} | {rating} | {safe} | {row['doom_score']:.3f} |")
    lines.append('')
    (OUTPUT_DIR / 'table5_doom_loop_intensity.md').write_text('\n'.join(lines))


# ══════════════════════════════════════════════════════════════════════════════
# PHASE 7: FORWARD PROJECTIONS — DOOM LOOP UNDER DEMOGRAPHIC SCENARIOS
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n{'=' * 80}")
print("PHASE 7: FORWARD PROJECTIONS")
print("=" * 80)

# Use UN OADR projections to project future doom loop intensity
# Get projection data from full panel
proj = full[['iso3', 'year', 'old_dep', 'gdp_pc_ppp']].copy()
proj = proj[proj['year'].between(2025, 2055)]
proj_countries = result['iso3'].unique() if len(result) > 0 else []

# Get latest fiscal data for each country (assume fiscal stance persists)
latest_fiscal = df.groupby('iso3').apply(
    lambda g: g.dropna(subset=['govt_debt_gdp', 'exp_rev_gap']).iloc[-1] if len(g.dropna(subset=['govt_debt_gdp', 'exp_rev_gap'])) > 0 else None
)
latest_fiscal = latest_fiscal.dropna(how='all').reset_index(drop=True)

proj_results = []
for iso in proj_countries:
    country_proj = proj[proj['iso3'] == iso].sort_values('year')
    fiscal = latest_fiscal[latest_fiscal['iso3'] == iso]
    if len(country_proj) == 0 or len(fiscal) == 0:
        continue

    debt = fiscal['govt_debt_gdp'].values[0]
    gap = fiscal['exp_rev_gap'].values[0]
    pension = fiscal['pension_spending_gdp'].values[0] if pd.notna(fiscal['pension_spending_gdp'].values[0]) else None

    for _, row in country_proj.iterrows():
        if pd.notna(row['old_dep']):
            # Simple doom score projection (static fiscal + projected demographics)
            oadr_norm = (row['old_dep'] - 3) / (40 - 3)  # normalize to 0-1 on plausible range
            debt_norm = (debt - 0) / (250 - 0)
            gap_norm = (gap - (-15)) / (30 - (-15))
            score = 0.35 * oadr_norm + 0.35 * debt_norm + 0.30 * gap_norm
            proj_results.append({
                'iso3': iso, 'year': int(row['year']),
                'oadr': row['old_dep'], 'doom_score': score
            })

proj_df = pd.DataFrame(proj_results)
if len(proj_df) > 0:
    # Show projection for key years
    for yr in [2030, 2040, 2050]:
        snapshot = proj_df[proj_df['year'] == yr].sort_values('doom_score', ascending=False)
        if len(snapshot) > 0:
            print(f"\n  Doom Score Projections — {yr}:")
            for _, row in snapshot.head(10).iterrows():
                print(f"    {row['iso3']:5s}  OADR={row['oadr']:.1f}%  Score={row['doom_score']:.3f}")

    # Write projection table
    lines = ['# Table 6: Doom Loop Projections (2030-2050)\n']
    lines.append('Projected doom loop intensity using UN demographic projections + current fiscal stance.\n')
    lines.append('| Country | 2024 Score | 2030 OADR | 2030 Score | 2040 OADR | 2040 Score | 2050 OADR | 2050 Score |')
    lines.append('|:---|---:|---:|---:|---:|---:|---:|---:|')

    for iso in proj_countries:
        current = result[result['iso3'] == iso]
        if len(current) == 0:
            continue
        current_score = current['doom_score'].values[0]
        row_parts = [f"| {iso} | {current_score:.3f}"]
        for yr in [2030, 2040, 2050]:
            snap = proj_df[(proj_df['iso3'] == iso) & (proj_df['year'] == yr)]
            if len(snap) > 0:
                row_parts.append(f" {snap['oadr'].values[0]:.1f} | {snap['doom_score'].values[0]:.3f}")
            else:
                row_parts.append(" — | —")
        lines.append(' |'.join(row_parts) + ' |')
    lines.append('')
    (OUTPUT_DIR / 'table6_doom_loop_projections.md').write_text('\n'.join(lines))


print(f"\n{'=' * 80}")
print("ALL PHASES COMPLETE")
print(f"Output tables written to: {OUTPUT_DIR}")
print("=" * 80)
